Intrinsic optical bistability of thin films of linear molecular aggregates: The 

one-exciton approximation 
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We perform a theoretical study of the nonlinear optical response of an ultrathin film consisting 
of oriented linear aggregates. A single aggregate is described by a Frenkel exciton Hamiltonian 
with uncorrelated on-site disorder. The exciton wavefunctions and energies are found exactly by 
numerically diagonalizing the Hamiltonian. The principal restriction we impose is that only the 
optical transitions between the ground state and optically dominant states of the one-exciton mani- 
fold are considered, whereas transitions to other states, including those of higher exciton manifolds, 
are neglected. The optical dynamics of the system is treated within the framework of truncated 
optical Maxwell-Bloch equations in which the electric polarization is calculated by using a joint 
distribution of the transition frequency and the transition dipole moment of the optically dominant 
states. This function contains all the statistical information about these two quantities that govern 
the optical response, and is obtained numerically by sampling many disorder realizations. We derive 
a steady-state equation that establishes a relationship between the output and input intensities of 
the electric field and show that within a certain range of the parameter space this equation exhibits 
a three-valued solution for the output field. A time-domain analysis is employed to investigate the 
stability of different branches of the three-valued solutions and to get insight into switching times. 
We discuss the possibility to experimentally verify the bistable behavior. 

PACS numbers: 42.65.Pc, 71.35.Aa; 78.66.-w 



I. INTRODUCTION 

Optical circuits make use of light to process informa- 
tion. They operate at the speed of light with almost no 
energy dissipation, unlike electronic analogs. Optical fi- 
bresi'^ and photonic crystal fibers^ have already found 
important applications in optical communications and 
optoelectronic devices. Implementing ultrafast optical 
sources and all-optical switches based on novel (quantum- 
confined) materials, such as organic thin films and quan- 
tum dotS)^ as well as silicon-based structureSf^ is now 
in progress. The realizability of a single-photon optical 
switch based on warm rubidium vapor has recently been 
demonstratedi^ 

A key element of any optical logic device is the opti- 
cal switch, which either passes or reflects the incoming 
light, depending on its intensity. One possibility to de- 
sign an optical switch is to utilize the phenomenon of op- 
tical bistability. Since the theoretical prediction of this 
effect by McCalF and its experimental demonstration by 
Gibbs, McCall, and Venkatesan* for a cavity filled with 
potassium atoms, an extensive literature, both theoret- 
ical and experimental, has developed on this topic (see 
Refs. ^9T0,11 for historical overviews and Ref. [13 for re- 
cent developments on optical instability in wide aperture 
laser systems). A generic optical bistable element ex- 
hibits two stable stationary transmission states for the 
same input intensity, a property which in principle opens 
the door to applications such as all-optical switches, op- 
tical transistors, and optical memories. 

Nonlinearity and feedback are the two necessary ingre- 
dients in order to enable optical bistable response of an 



optical system. The former can be provided, e.g., by a 
saturable medium, while a cavity (mirrors) can serve to 
build up a feedback. This arrangement has been used 
in the first demonstration of controlling light with lights 
Sometimes, however, the nonlinearity itself plays the role 
of the feedback. Here, bistability is an intrinsic property 
of the material; no external feedback, like a cavity, is 
needed. Thus, mirrorless (or cavityless) optical bista- 
bility can be realized, which is even more advantageous 
from the viewpoint of designing all-optical devices. Dur- 
ing the past decade, this type of bistability has been ob- 
served in a variety of inorganic materials heavily doped 
with rare-earth ions i^^i^^'^^i^° In Refs. [H, a population 
dependent dipole-dipole interaction in ion pairs has been 
put forward as a nonlinearity and feedback mechanism to 
explain the effect. This interpretation has been debated 
in a number of paperS ]^'*'^^i^^'^^ >^^'^^'^" 

Another class of materials, promising from the view- 
point of all-optical manipulation of light, are molecu- 
lar aggregates and conjugated polymers. These systems 
commonly exhibit narrow absorption bands and suppres- 
sion of exciton-phonon coupling, superradiance and giant 
optical nonlinearities, fast collective optical response and 
efficient energy or c harge t ransport (see for an overview 
Refs. IHiiEllii^IllSi), which are ingredients nec- 
essary to design optoelectronics or all-optical devices. 
Molecular aggregates and conjugated polymers have al- 
ready been used to fabricate light emitting diodes^i and 
organic solid-state lasers.-- 

One particularly interesting effect, which has already 
received a considerable amount of theoretical discussion, 
but still awaits experimental realization, is the mirrorless 
optical bistability of a single molecular aggregate2£ or 
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an assembly of molecular aggregates] i'^^' The bistable 
behavior of a single linear aggregate consists of a sud- 
den switching of the aggregate's excited state population 
from a low level to a higher one upon a small change of the 
input intensity around a critical point. The effect orig- 
inates from a dynamic resonance frequency shift, which 
depends on the number of excited monomers in the ag- 
gregate. The origin of this shift lies in the quasi-fermionic 
nature of Frenkel excitons in one 

dimension,Mi25i36 r^j^^^ 

nonlinearity plays the role of intrinsic feedback, neces- 
sary for bistability to occur. There exists, however, a re- 
striction on the aggregate length: an aggregate exhibits 
bistable behavior only if its coherence length is larger 
than the emission wavelength, which makes experimen- 
tal realization problematic. 

An assembly of molecular aggregates arranged in an ul- 
trathin film geometry (with the film thickness small com- 
pared to the emission wavelength) may display intrinsic 
optical bistability governed by another mechanism, where 
the density of molecules becomes the driving parameter. 
The same mechanism holds for an ultrathin film of homo- 
geneously broadened two-level systems. "^^ When the den- 
sity in the film is high enough, the on-resonance refractive 
index can get sufficiently large to totally reflect an incom- 
ing field of low intensity. Then the incoming field is al- 
most completely compensated by a secondary field of op- 
posite phase, which is generated by the aggregate dipoles. 
The dipole-induced field is bounded in magnitude, mean- 
ing that this picture only holds if the incoming field in- 
tensity is smaller than a certain value, determined by 
the density of aggregates. When this value is exceeded, 
the aggregates become saturated, which suppresses the 
dipole-induced field and abruptly changes the (nonlin- 
ear) refractive index and transmittivity of the film. The 
field produced by the aggregate dipoles plays the role of 
intrinsic feedback. The output field depends nonlinearly 
on the input field of the film. 

In Refs.lsil-lssia thin film arrangement of oriented lin- 
ear J-aggregates was considered, where the localization 
segments of a single disordered aggregate were modeled 
as independent homogeneous chains of fluctuating size. 
Each segment was considered as a few-level system, with 
an individual ground state and one or two excited states 
corresponding to the dominant optical states of the seg- 
ment. Within this framework, both the ground state to 
one-exciton2i and one-to-two-— exciton transitions were 
taken into account, and bistable behavior was found in a 
certain region in the parameter space. 

The approach used in Refs. 'Sl'andlsi assumed full cor- 
relation of fluctuations of the lowest exciton energy and 
the transition dipole moment, taking both magnitudes as 
solely depending on the segment size. The real picture, 
however, is quite different i^i^S In practise the optical re- 
sponse of J-aggregates is strongly affected by disorder in 
the molecular transition energies. The band-edge of the 
exciton energy spectrum of such a disordered aggregate 
is formed by states that are localized on segments with 
small overlap. The lowest state of a segment is optically 



dominant, whereas the other states have a much smaller 
oscillator strength. The energy of the lowest state is not 
correlated with the size of the segment; it is determined 
by uncorrelated well-like fluctuations of the site poten- 
tial4^ Therefore, the optically dominant states of non- 
overlapping segments can be arbitrarily close in energy, 
having at the same time completely different transition 
dipoles.*^ In other words, the transition dipoles and en- 
ergies of the relevant states turn out to be uncorrelated 
rather than correlated. 

In this paper, we exploit the two-level model, imple- 
mented in Refs. [3l| andls^. to describe the film's optical 
response. However, unlike Refs. [U and [H, we will ac- 
count properly for the statistical fluctuations of the tran- 
sition dipole moment and the transition energy, as they 
appear after diagonalizing the Frenkel exciton Hamilto- 
nian with uncorrelated on-site disorder. We calculate the 
joint probability distribution of these quantities and use 
it to compute the electric polarization of the film, which 
features in the Maxwell equation for the field. The ag- 
gregate segment dynamics is described within the 2x2- 
density matrix formalism. We derive a novel steady-state 
equation for the output field intensity as a function of the 
input intensity in terms of the joint probability distribu- 
tion of the energy and the transition dipole moment. On 
this basis, the bistability phase diagram of the film is 
calculated. The critical parameter for bistability to oc- 
cur turns out to be different (larger) than that found in 
Refs. Isi]. By numerically solving the truncated Maxwell- 
Bloch equations in the time domain, we study the stabil- 
ity of the different branches of the three- valued solution 
for the output field intensity. The calculation of an opti- 
cal hysteresis loop (an adiabatic up-and-down-scan of the 
field) demonstrates that only two of them are stable. A 
new element in the paper is that we also analyze switch- 
ing time between both stable branches, and show that it 
slows down dramatically close to the switching point. 

The outline of this paper is as follows. In section 
we present the model and mathematical formalism. Sec- 
tion |TTT] deals with the linear regime of the transmission. 
The steady state equation for the output intensity in 
the nonlinear regime is derived in Sec. IIVI In Sec. |Vl 
the stability of different branches is considered, together 
with a study of the switching time. In Sec. IVII we dis- 
cuss the possibility to achieve optical bistability using 
J-aggregates of polymethine dyes. Section IVIII summa- 
rizes the paper. Finally, in the Appendix we address the 
effect of interference of the ground state to one-exciton 
transitions, originating from the fact that excitons are 
born from the same ground state, with all monomers be- 
ing unexcited. 



II. MODEL AND FORMALISM 

We aim to study the transmittivity of an assembly 
of linear J-aggregates arranged in a thin film geometry 
(with the film thickness L small compared to the emis- 
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sion wavelength A' inside the film). All aggregates are 
aligned in one direction, parallel to the film plane. Such 
an arrangement can be achieved, e.g., by spin-coating.— 
The limit of L ^ A' allows one to neglect the inhomo- 
geneity of the field inside the film. The aggregates in 
the film are assumed to be decoupled from each other. 
This finds its justification in the strongly anisotropic na- 
ture of the system we have in mind. As we will see later 
(Sec. VI), films of interest for bistability should have a 
molecular density of the order of 10^^ cm~'^. With a typi- 
cal separation of 1 nm between molecules within a single 
aggregate, this implies that neighboring aggregates are 
separated by 10 nm. Thus, the dominant dipole-dipole 
interactions between molecules of different chains are a 
factor of 1000 weaker than those within chains. As a 
consequence, we expect that the former interactions will 
merely result in small shifts of resonance energies, away 
from the single-chain exciton energies considered below. 

On the other hand, the effect of interactions of the ag- 
gregate molecules with the surrounding host molecules 
is important, because as a consequence of the usually 
inhomogeneous nature of the host media, they lead to 
disorder in the molecular transition energies and in the 
molecular transfer integrals, both of which give rise to 
localization of the exciton states on segments of the ag- 
gregates. Finally, thermal fluctuations in the environ- 
ment result in intraband scattering of the excitons that 
causes two effects: equilibration of the exciton popula- 
tion and homogeneous broadening of the exciton levels. 
In this paper, we neglect the former effect. This finds its 
justification in many experimental studies, which have 
shown that the fiuorescence Stokes shift of J-aggregates 
of cyanine dyes usually is very small 4^'^'^'^) 

A. A single aggregate 

We model a single aggregate as a linear array of N two- 
level monomers with parallel transition dipoles. In this 
paper, we restrict ourselves to optical transitions between 
the ground state an the one-exciton manifold, described 
by the Frenkel exciton Hamiltonian 

N N 

Ho = ^ e„|n)(n| -f ^ J„„ \n){m\ , (1) 

n— 1 n.m 

where \n) denotes the state in which the nth site is ex- 
cited and all the other sites are in the ground state and e„ 
is the excitation energy of site n. The e„ are taken at ran- 
dom and uncorrelated from each other from a Gaussian 
distribution with mean eg (the excitation energy of an iso- 
lated monomer) and standard deviation a. The transfer 
interactions Jnm are considered to be of dipolar origin 
and non fiuctuating: Jnm = ~J/W ~ '^P {Jnn = 0). 
Here the parameter J represents the nearest-neighbor 
transfer interaction, which will be chosen positive (as is 
appropriate for J-aggregates). The exciton energies 
{v = 1,...,N) and wavefunctions = "^n^ifunln) , 



are obtained as eigenvalues and eigenvectors of the N xN 
Hamilton matrix Hnm = {n\H\m). 

From the set of exciton states we only take into 
account the optically dominant states which, for J > 0, 
reside in the neighborhood of the low-energy bare band 
edge at eq = eq — 2.404J. These states are located 
at different segments of the aggregate, which overlap 
weakly, and have a wavefunction with no node. There- 
fore, they are called s-like states. To find all such states, 
we use the selection rule proposed in Ref. [H. It reads 
> ^0, where we set Co = 0.8. This rule 
selects states with a wavefunction consisting of mainly 
one peak. From now on, the state index v will count only 
such s-like states. The number of these states is roughly 
equal to N/N*, where TV* is their typical localization size. 
We assume that the vibration-induced coherence length 
of excitons is much larger than the disorder-induced lo- 
calization length, a condition that can be fulfilled at low 
temperature^^ In this limit, the exciton eigenstates 
form a good basis. 

The above picture implies that an aggregate is mod- 
eled as a set of independent segments, each of which has 
its own ground state |0) and an s-like excited state \v). 
The optical transition between these states is governed 
by the segment dipole operator d„ = do(|0)(j^| -I- |j^)(0|), 
where cIq is the transition dipole moment of a monomer. 
The corresponding transition dipole moment of a seg- 
ment is calculated as = do'^^ipun = rfoM^, where 
fii, — if^n is the dimensionless transition dipole mo- 
ment. 

The optical dynamics of a segment is described in 
terms of the 2 x 2-density matrix (p^,^ , p^Q , p* q , poo ) which 
obeys the Bloch-like equations (see the Appendix) 

Puu = -luPuv + idu^ {pou - Puo) , (2a) 
Puo = - ii£u + r^) p^o - iduS{puv - Poo) , (2b) 

Poo + p^i. = 1 . (2c) 

Here we set the Plank constant h — 1 and introduced 
the following notations: 7^/ = 7o|pi^P is the radiative 
rate of the exciton state v (70 being the monomer radia- 
tive rate), and F^ = ^7^ -I- 71^0 is the dephasing rate of 
the state v, which includes a pure dephasing term, 71,0- 
Finally, E is the total electric field inside the film (see be- 
low). Owing to the disorder, the transition energy e^, the 
relaxation constant F^, and the transition dipole moment 
Pi, are stochastic variables, which differ from segment to 
segment. Because of the fiuctuations in e^, F^, and d,j, 
the density matrix elements pi^i,, p^o, and poo fluctuate 
as well. 



B. The Maxwell equation 

In this section, we specify the field £ which enters 
Eqs. It consists of two contributions: the incoming 
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field Si and a part produced by the aggregate dipoles. 
The incoming field is considered to be a plane wave 
£i ~ Ei{x,t) cos{kiX ~ LOit) with a frequency LOi and an 
amplitude Ei{x, t), normally incident and polarized along 
the aggregate transition dipoles. Under these conditions, 
all the vectorial variables (transition dipole moments, in- 
coming and outgoing fields, and the field inside the film) 
can be considered as scalars. The amplitude Ei{x,t) is 
assumed to vary slowly on the scale of the optical period 
2TT/uJi and wavelength = 2i:/ki. 

We assume without loss of generality that the film is 
located in the ZY plane {x = 0). Then the total field at 
a; = (inside the film) is given by^i^ 



. . 27rL . 
£ = £^ V , 



(3) 



where P is the electric polarization of the film (electric 
dipole moment per unit volume), the dot denotes the 
time derivative, and c stands for speed of light. The sec- 
ond term in the right hand side of Eq. ([3]) represents the 
field produced by the dipoles in the film, emitted perpen- 
dicular to the film in both directions. The part propa- 
gating to the left is the reflected (plane wave) field, given 
at x = by f r = — (27rL/c)P, while the part propagating 
to the right is the emitted (also plane wave) field, which 
forms, together with the incident field £i , the transmitted 
signal, determined at cc = by Eq. ([3]). 

The electric polarization V is calculated as follows. 
First, we introduce the expectation value of the dipole 
operator of an aggregate, d = d^Y^^^s Mi^CPi-o + Pou), 
where the summation is performed only over the s-like 
states of the aggregate. Furthermore, this value is aver- 
aged over a physical volume V , containing M aggregates, 
which, in fact, is equivalent to obtaining the average {d) 
over disorder realizations. After that, V is obtained by 
multiplying (d), with the number density M/V of the ag- 
gregates. The final formula for the electric polarization 
reads: 



V = dono^ 



(4) 



Here, tiq = NM/V is the number density of monomers. 



Ns ~ ^ X^i/es ^1 ^ normalization constant (having the 
meaning of the average number of s-like states in an ag- 
gregate), and pio{£,^,t) is the off-diagonal density ma- 
trix element, where the indices and 1 label the ground 
and the excited s-state of the segment, respectively. In 
our present formulation this element, as well as poo and 
pii, are ordinary (not stochastic) functions of e and fi; 
which formally follow from solving Eqs.®. All stochastic 
aspects of the segment's properties are taken into account 
through the function Gsis, fi), which represents the joint 
probability distribution of the transition energy e and the 
dimensionless transition dipole moment fi of the segment. 
The latter is defined as 
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FIG. 1: (a) The joint probability distribution Qa{e,p) of the 
transition energy e and dimensionless transition dipole mo- 
ment /i for s-like states on localization segements, obtained 
for a disorder strength a = O.IJ according to Eq. ((5)1. We 
used chains of length N = 500 with the monomer transi- 
tion energy eo = 0. The sampling was performed over 300 
000 disorder realizations. Contour lines correspond to 10% 
of the peak value of the distribution, (b) - The absorption 
spectrum Aa{s) ~ J dp p?Qs{e, p). (c) - The distribution 
Ms{p) ~ J deQs{e,iJ.) of the transition dipole moment p. 
The solid lines represent the results of calculations, whereas 
the open circles are fits by a Gaussian. 



It is worth to notice that at a given disorder strength 
a, Ns scales linearly with the aggregate size N. Hence, 
the ratio Ng/N in Eq. ^ is iV- independent. From our 
simulations we found that N^/N = 0.074(cr/ J)°-^, which 
nicely agrees with the disorder scaling of the typical lo- 
calization size iV*j^ 

After the fJs-distribution is obtained by straightfor- 
ward sampling of a sufficient number of disorder realiza- 
tions, one can easily calculate the two important quanti- 
ties: As{e)^N-\J2^^^filS{e~e,)) ^ Jdfifi^gs{e,ti), 
which represents the absorption spectrum, not account- 
ing for homogeneous broadening (i.e., close to the zero- 
temperature spectrum), and Ais{p) = X^i/gs 
fJ.^)) = J (ie^s(e,/i), which represents the probability 
density of the transition dipole moment. As we are 
mostly interested in the limit of dominating ingomoge- 
neous broadening, we will refer from now on to As{s) as 
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to the absorption spectrum, assuming that its half width 
at half maximum (HWHM) a* is larger than the homoge- 
nous HWHM (resulting from Yy). 

An example of the distributions Qs(e^[i)^As{s), and 
A^s(/i) computed for an ensemble of chains with N ~ 500 
and a disorder strength g = 0.1 J, is depicted in Fig. [1] 
[panels (a),(b), and (c), respectively]. Because Qs{e,fi) — 
Gs{£, — Ai), only /x > is considered in the plots. 

Note that in our model, the absorption spectrum As (e) 
is almost symmetric with respect to the peak position, 
except the tails, which show a small asymmetry. It can 
be fitted well by a Gaussian, unlike the case when all the 
exciton states are taken into account. The latter gives 
rise to a Lorentzian high-energy tail of As{s), reproduc- 
ing the asymmetric lineshape commonly seen in exper- 
iments. The shape of the A^s-distribution can also be 
fitted by a Gaussian, but with a lesser accuracy than the 
absorption spectrum. The distribution Qs{£,^) exhibits 
interesting scaling properties with regard to the disorder 
strength a. A detailed study will be presented elsewhere. 



C. Truncated Maxwell-Bloch equations 

To proceed we seek the solution of Eqs. ^ in the 
standard manner: we set pio — ~{i/2)Rexp (—iujit) and 
£ = (l/2)i?exp (—iujit) + c.c, where the complex ampli- 
tudes R and E vary slowly on the time scale iir/uji, and 
we use the rotating wave approximation. It is straight- 
forward to arrive at a set of truncated equations for the 
populations pn of the one-exciton states, and the ampli- 
tudes R of the off-diagonal density matrix elements, and 
the field f2 = d^E (in frequency units): 

Pii = -7P11 - i^R* + ^*R) , (6a) 



i? - - - Ao) + r] i? + pn{pn ~ Poo) , (6b) 



N f 

Vr^ I dAdfig{A,p)fiR 



(6c) 



where A — Aq — e — iOt is the frequency detuning between 
the exciton transition and the incoming field, which is de- 
composed into two parts: A = e — Sq and Aq = uJi — Sq 
indicating, respectively, the frequency detuning of the ex- 
citon state and the incoming field from the exciton band- 
edge frequency Sq = eo — 2.404J. 

The constant Tn = 2'Knodo^kL is an important pa- 
rameter of the modeli^i^^i^ The physical meaning of 
Til can be explored by rewriting it in the form Tn = 
^jonoL{X/2)^ , where 70 = 4dQa;'^/(3c'^) is the monomer 
spontaneous emission rate and hqL is the surface density 
of monomers. The quantity noL(A/2)^ can be interpreted 
as the number of monomers in a (A/2)^-square that oscil- 
late in phase. In other words, Tn can be considered as the 
radiative rate of a single monomer, 70, enhanced by the 



number of monomers within a (A/2)^-square4^ Tr gov- 
erns the Dicke superradiance of a thin filmj^i^ as well 
as the collective radiative damping in the linear regime 
(see the next section). Therefore it is usually referred to 
as the superradiant constant. 

The set of equations (O together with the normaliza- 
tion condition (Pc)l and the definition ^ form the basis 
of our analysis. In the remainder of this paper, we will be 
particularly interested in the dependence of the transmit- 
ted field intensity |r2|^ on the input field intensity jilip 
following from these equations. 



III. LINEAR REGIME 

In order to get insight into the effect and interplay 
of the parameters that govern the bistability, we first 
consider the linear regime of the system. We assume that 
a weak input field Cli = const is switched on at i = 0, 
weakness implying that the depletion of the ground state 
population can be neglected. Thus, we set poo{t) = 1 and 
Piiit) = 0, which linearizes Eqs. ([6]), 



R = - [i{A - Ao) + r]R- fin , 



Ns 



(7a) 



n = n, + TR^ J dAdngsiA,p)pR. (7b) 

These equations can be solved easily in the Laplace do- 
main. The solution for the Laplace transform of the 
transmitted field 51 reads 



l + Tfl^ / dAdpgs{A,p)p^ 



>+[z(A- Ao) + r] 



(8) 



where p denotes the Laplace parameter. To evaluate this 
expression, we neglect the p-dependence of T. Then the 
integral over p, gives, by definition, the absorption spec- 
trum As{A). The latter is normalized now to Fg/Ng, 
where Fg — ( X^^igs f^i) ^^'^ average total oscillator 
strength of the s-like states per aggregate. To perform 
the integration over A explicitly, we replace ^^(A) by a 
Lorentzian centered at A* and with a width a*: 



MA) = ^ ^ 



1 



Ng IT [(A- A*)2+(T*2] 



(9) 



(in all our numerical results, we do not invoke this 
approximation and keep the exact form of the Gg- 
distribution, i.e., of the absorption spectrum). With this 
substitution, the result of the integration over A reads: 



p + i{A* Ao) + r + cr* + ri 



(10) 



where we introduced the renormalized superradiant con- 
stant Tn = {Fg/N)rii. As the total oscillator strength 
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of s-like states Fg < N, the ratio Fs/N < 1. We also 
note that F + tr* denotes the total (homogeneous plus 
inhomogeneous) dephasing rate. 

Finally, by assuming = const, which corresponds to 
fli — ^i/s in the Laplace domain, we obtain the following 
time-domain behavior of the transmitted field 



n = 



i(A*-Ao)+F + o-* 
i(A* - Aq) + F + cr* + 



i{A* - Ao) + r + (7* + Tb. 
X exp -i{A* - Ao)t- (T + a* +rii)t 



(11) 



As is seen from this equation, the field in the film, Q, 
reaches its steady-state value (given by the first term in 
the right-hand side) after a time l/(r + cr* + F/j). If the 
dephasing dominates the relaxation of the dipolcs, i.e., 
if F + cr* ^ Tji, the steady state limit of the opposing 
dipole field, given by -^2,fii/[i(A* - Ao) + F + a* + Ffl], 
is small in magnitude compared to the incoming field fli. 
As a consequence, the field inside the film « fJj. In 
this limit, one finds a high film transmittivity. 

When Tu ^ F -I- ct* the superradiant damping drives 
the relaxation. Now the film dipoles, having sufficient 
time to respond collectively, can produce an opposing 
field -rjjf fl/|i(A* - Ao) + T + a* + TrI of magnitude 
« 51 j. This field almost totally compensates the incoming 
field, and results in a low magnitude of the field inside 
the film, \n\ - a,\i{A* - Aq) + T + a*\/rR < fli, and, 
consequently, in a low film transmittivity. Switching to a 
high transmission state now requires a field intensity Qi 
that saturates the system. In this case we can see optical 
bistable switching (see the next section). 

From the above, it is clear that the interplay of su- 
perradiance and dephasing determines the linear trans- 
mittivity of the film. Hence, the ratio Tji/{r -I- cr*) is 
an important parameter of the model. In the theory of 
bistability it is often referred to as the cooperative num- 



IV. STEADY-STATE ANALYSIS 

A. Bistability equation 

In this section, we consider the steady-state regime, 
when we set Vti = const and R = pu = 0. It is a matter 
of simple algebra to derive the following equation for the 
output intensity |ilp: 



l + Tfl ^ I dAdtiti^g,iA,^i) 



F-z(A- Ao) 



(A-Ao)2+F2 + |fl|2F/7o 



(12) 



Formally, Eq. p2)) differs from the one found previ- 
ouslySi by the small factor Ng/N. This smallness, how- 
ever, is compensated by the A^^-scaling of the integral 
in the latter is proportional to Fg/Ng ^ 1 (see the 
preceding section). Thus, the actual numerical factor in 
Eq. (fT2|) is on the order of Fg/N. Numerically, we found 
that Fg/N depends only weakly on the disorder strength 
cr, lying within an interval 0.75 < Fg < 0.83 when the 
disorder strength a ranges from to 0.5J. This means 
that the linear optical response in a system with static 
disorder is dominated by the s-like states, independent 
of the disorder. We stress that, unlike previous worksf^ 
Eq. (|12p properly accounts for the joint statistics of the 
transition energy and the transition dipole moment via 
the ^/^.-distribution. 



B. Phase diagram 

Numerical analysis shows that Eq. can have three 
real roots in a certain region of the parameter space 
{Tr, cr*). In other words, our model can exhibit bistable 
behavior. In all simulations, we used linear chains of 
N = 500 sites and 70 = 2 x 10~^J (appropriate for 
monomers of polimethine dyes). The dephasing con- 
stant was considered not fluctuating^ and was set 
to 7^0 = 5OO70. 

Several examples of the S'-shaped input-output charac- 
teristics calculated for the disorder degree cr — O.IJ are 
shown in Fig. [2] for an input field that is resonant with the 
absorption maximum. We use the dimensionless intensi- 
ties /in = |fljp/(7ocr*) and /out = |flp/(7ocr*), which is 
convenient because the HWHM of the absorption spec- 
trum cr* is an experimentally measurable quantity. Panel 

(a) shows how the input-output characteristics change 
when Tr is below, at, or above its critical value. Panel 

(b) shows the input-output characteristics when the field 
is tuned through the resonance. 

At a given disorder strength cr, the minimal value of 
the superradiant constant Tr needed for optical bista- 
bility (the critical value F^) depends on the detuning 
Aq. Fig ure [Slja) explicitly demonstrates this effect: F^ 
is almost constant within the absorption band, whereas 
it clearly increases outside it. Panel (b) shows the Aq- 
dependence of the critical switching intensity /^ of the 
incoming field at the bistability threshold. This is the 
lowest intensity at which the film can switch, when the 
field is tuned at Aq, and when the superradiance constant 
Tr = F^(Ao). The data presented here is obtained for 
the disorder strength cr = O.IJ. 

As is seen from Fig.[3Ua), there exists a detuning Aq^*, 
referred to as the optimal one, at which F^ takes its mini- 
mal value. The detuning is optimal if the imaginary term 
in Eq. vanishes: this term opposes a three- valued so- 
lution for the output field. For a symmetric absorption 
band, the optimal detuning corresponds to the incoming 
field being resonant with the absorption maximum. In 
our case, owing to a small asymmetry of the absorption 
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FIG. 2: Examples of the input-output characteristics, demon- 
strating the occurrence of three- valued solutions to Eq. (I12|) . 
In simulations, chains of A'' = 500 sites and a disorder 
strength a = O.IJ were used, corresponding to a HWHM 
a* — 0.0156J. (a) - The results obtained for different super- 
radiant constants Tr at the optimal detuning Ag''* = —2.42 J, 
which corresponds to an incoming field which is resonant 
with the absorption maximum. The open circles, dotted, 
and solid curves represent, respectively, the data calculated 
for r_R — I6.6I0"* (the bistability threshold for a = O.IJ), 
r_R = 11.52a* (below the bistability threshold), and 27.12(t* 
(above the bistability threshold), (b) - The results obtained 
for Tn — 16.61(7* and various detunings Aq. The dotted and 
solid curves represent, respectively, the data calculated for 
Ao = Aq''' — a", and Aq''' + a*. The open circles show the 
same data as in panel (a). 



FIG. 3; (a) - Dependence of the critical superradiant constant 
on the detuning Aq (solid line) calculated for the disor- 
der strength a = O.IJ. The dashed line shows the absorption 
spectrum (absorption only due to s-states). The dotted hor- 
izontal line indicates calculated for the optimal detuning 
^opt _ „2 42 J. (b) - Dependence of the switching inten- 
sity Jin on the detuning Aq calculated at the corresponding 
bistability threshold, i.e., with F^ given in the panel (a). 

tion is the same as for a thin film of homogeneously 
broadened two-level systems, only with the renormalized 
cooperative number Tr/T* — (Tr/T*){Fs/N), where 
Fs/N = 0.83. Bearing in mind that the critical value of 
the ratio Tr/T* is equal to 8,^1 we recover F^/F* = 9.64. 



band [see Fig.Wl.b)], Aq^* = -2.42J is shifted slightly to 
the blue from the position of the absorption peak. 

We calculated F^ as a function of the HWHM a* for 
the optimal detuning. The result is shown in Fig.[?l The 
plot represents, in fact, the phase diagram of the optical 
response: below the curve, the output-input character- 
istic of the film is always single- valued (stable), while 
- depending on the detuning - it can become three- 
valued (bistable) above it. The nonmonotonic behav- 
ior of F^ at small magnitudes of a*, presented in the 
panel a, is simply explained by the fact that the disorder- 
induced (inhomogeneous) broadening becomes smaller 
than the homogeneous one a* < F*, where F* is defined 
as F* = J diideTQs{e, n). The ratio Tr/T* is now the 
relevant parameter, governing the occurrence of bista- 
bility. The panel (6) shows the (7*-dependence of F^ 
replotted in units of F*, which is monotonic. When 
cr ^ 0, the ratio F^j/F* ^ 9.64. This value is de- 
duced from Eq. Indeed, in the limit of cr ^ 
we can move the Lorenztian outside the integral and 
use / dAdfj, fi^QsiA, ^) = Fg/Ng. The resulting equa- 



C. Spectral distribution of the exciton population 

More insight into what occurs at the switching thresh- 
old is obtained by studying the population distribution 

rn(A) = J d/.g,(A,/i)pn(A,/i) , (13) 

with pii the steady-state solution of Eqs. ([6]). This distri- 
bution enables us to visualize the relevant spectral range 
around the switching point. 

In Fig. [5l we plotted rii(A) calculated for the opti- 
mal detuning Aq''' and Tr = 27.12cr* (above the criti- 
cal value F^). Panels (a) and (b) show the results ob- 
tained for the incoming field intensities Lm = ^i/ilo<^*) 
below and above the switching threshold, respectively. 
Below the switching threshold, only a relatively narrow 
spectral region around Aq''' acquires population. This is 
because, in spite of the intensities of the incoming field 
/in = 3.33,64.34, and 81.78 being far above the satu- 
ration value, the intensity of the field inside the film. 
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FIG. 4: (a) - Phase diagram of the bistable optical response of 
a thin film in the (Tr, a*)-space obtained by solving Eq. (|12|) 
for A = Ag^'. The open circles represent the numerical data 
points, whereas the solid line is a guide to the eye. Above 
(below) the solid line the film behaves in a bistable (stable) 
fashion. The solid line itself represents the a"*-dependence 
of the critical superradiant constant F^, calculated for the 
optimal detuning Aq''*, i.e., when the incoming field is tuned 
to the absorption band maximum. This gives the minimal F^ 
for each a* . (b) - The same data points as in the panel (a), 
only replotted as a function of F* , where F* is the mean value 
of the relaxation constant F. 




FIG. 5: Population-distributions rii(A) (solid curves), calcu- 
lated according to Eq. (Hg for cr = O.IJ and F^ = 27.12cr*, 
with the optimal detuning Aq^' = — 2.42J indicated by the 
vertical dashed line. Open circles show the absorption spec- 
trum As(A). Panel (a) represents rii(A) below the upper 
switching threshold. The plotted distributions were calcu- 
lated for the input intensities Jin = /{'joa*) = 3.33, 64.34, 
and 81.78 (from bottom to top). Panel (b) shows rii(A) 
above the upper switching threshold. In the inset, the depen- 
dence of the full width at half maximum (FWHM) of rii(A) 
on Jin is plotted in units of the FWHM of the absorption 
spectrum. 



ing the one-to-two exciton transitions. This work is now 
in progress. 



V. TIME-DOMAIN ANALYSIS 



A. Hysteresis loop 



/out = ^^1/(700"*) — 0.025, 0.5, and 1.5, is below or on 
the order of it. For these intensities, the one-exciton ap- 
proximation, with only one s-like excited state considered 
in each localization segment, is reasonable. 

Figure EJb) represents the population distribution 
rii(A) after switching, when the field inside the film 
/out exceeds the switching threshold and becomes much 
larger than the saturation magnitude. In this limit, 
we can replace pii{A, fi) in Eq. by 0.5 and get 

rii(A) = 0.5 / duGsiA^n) ^ 0.52?^ (A), where P^(A) is 
the density of s-like states. The latter is plotted in Fig. [5] 
(b) by the solid line and appears to be wider than the 
absorption band. For such field intensities, it is likely 
that the two-level model should be corrected by includ- 



It is well known that the S-shaped output-input depen- 
dence and, as a consequence, the existence of two switch- 
ing thresholds results in optical hysteresis. ^^'^^ To investi- 
gate this, we numerically integrated Eqs. ^ while slowly 
sweeping up-and-down the input intensity /in above the 
bistability threshold (Ffl > F'^). The result for the trans- 
mitted intensity /out is shown in Fig. [6] by the solid curve 
with arrows. The parameters used in the calculations 
are specified in the figure caption. The input field inten- 
sity was swept from zero to 110 and back to zero. The 
open circles indicate the steady-state solution obtained 
by solving Eq. for the same set of parameters. 

As can be seen from Fig. [51 the solid curve almost per- 
fectly follows the lower and upper branches of the steady- 
state three- valued solution, nicely demonstrating the op- 
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FIG. 6: An example of the stable optical hysteresis loop of 
the transmitted intensity lout ~ |rip/(7oo"*) (the solid curve 
with arrows) obtained by numerically solving Eqs. (|6]l for a 
linear sweeping up-and-down of the input field intensity /in ~ 
P/(7oo"*)- The sweeping time is 3000/(T* . The open circles 
represent the steady-state solution, Eq. (|12|) . The calculations 
were performed for the following set of parameters: r2 = 
5OO70, a = O.IJ, Tr = 27.12(7*, and Ao = AqP* = -2.42J. 



tical hysteresis. The intermediate branch is not revealed, 
which is clear evidence of its instability. Note also that 
switching from the lower branch to the upper one occurs 
for an input field intensity larger than the critical value. 
This indicates that when the input field intensity is only 
slightly above the switching intensity, the response of the 
film slows down. A much less pronounced but similar ef- 
fect can be observed at the lower switching threshold, 
where the field switches from the upper branch to the 
lower one. This is consistent with our study of the relax- 
ation time presented below. 



FIG. 7: Kinetics of the transmitted field intensity lout = 
jr2|'^/(7o(J*) approaching its stationary value (dashed line) af- 
ter the incident field with intensity Jin = |r2ip/(7ocr*) — 150 
is turned on abruptly ai t = 0. The value Jin ~ 150 ex- 
ceeds the upper switching threshold /;„ — 82.16. The other 
parameters were chosen as in Fig. [S] 
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B. Switching time 

Of great importance from a practical point of view, is 
the relaxation time r which is required for the output 
intensity to approach its steady-state value after the in- 
put intensity has changed. If this time is much shorter 
than the characteristic time of changing the input in- 
tensity, then the output signal will adiabatically follow 
it, remaining all the time close to the steady-state level. 
Only in the limit of short r, an abrupt switching from 
low to high transmittivity can be realized. This espe- 
cially concerns the region in the vicinity of the switching 
thresholds (see Fig. [5]). In other words, the relaxation 
time T limits the usage of the optical bistable element as 
an instantaneous switcher. 

Motivated by the above observations, we performed a 
study of the relaxation time t. Figure [7| shows an ex- 
ample of how the transmitted field intensity approaches 
its stationary value when an input field intensity with 
a value of 1^ = 150 is instantaneously switched on at 
t = 0. This field is above the upper switching threshold 



FIG. 8: Relaxation time r as a function of the excess input 
intensity /in — /in at the upper switching threshold (indicated 
by the vertical dotted line), r was calculated by turning on 
abruptly the incoming field a,t t — 0, and waiting until the 
transmitted field intensity /out approaches its steady-state 
value (for more details, see the text). The open circles show 
the numerical results, while the solid line represents a best 
power-law fit given by Eq. 1141 The calculations were per- 
formed for the set of parameters of Fig. (6) 



= 82.16. Calculations were carried out for the set 
of parameters of Fig. [6l As is observed, for the set of 
parameters used, the output intensity stays low during a 
waiting time of about 20/cr*, before it rapidly (on a time 
scale much shorter than 20/cr*) increases to its steady 
state value. This behavior allows one to define r as the 
time which the output intensity takes to reach its first 
peak (17.3/(7* in the current example). 

Using the above definition, we calculated the relax- 
ation time r as a function of the excess input intensity 
/in — at the upper switching threshold. The results 
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are plotted in Fig. [8] As is seen, r drastically increases 
when /in gets closer to The numerical data (open 
circles) is well fitted by the formula 

r = 870(/in-/f„)^°-'' , (14) 
shown by the solid curve. 

VI. DISCUSSION OF DRIVING PARAMETERS 

To get insight into the possibility to realize optical 
bistable behavior for a film of J-aggregates, we con- 
sider the typical parameters for this type of systems. 
First, we estimate the superradiant constant F/j — 
(3/87r)7onoA^L, considering the low-temperature exper- 
imental data of J-aggregates of polymethine dyes. For 
these species, typically, 70 ~ (1/3) ns^^ and A « 600 
nm i^'^i'^^i^^i^°i^^i^^ With this in mind and choosing L = 
\/2tt (or kL = l),we obtain the following estimate: 
Tr 10~^^no cm^ ps~^. This value for L is easily acces- 
sible with the spin-coating method^ and guarantees the 
applicability of the mean-field approach for the descrip- 
tion of the thin film optical response. '^^ The typical width 
of the J-band of polymethine dyes is on the order of sev- 
eral tens of cm~^ or approximately 1 ps~^ (in frequency 
units) Thus, for the set of parameters we 
chose, the number density of molecules ng must be on 
the order of 10^^ cm~'^ to get the ratio Tn/a* required 
for bistability to occur. This concentration is usually 
achieved in spin-coated films. 

Another option to adjust the parameters favoring 
bistability is to consider J-aggregates composed of 
monomers with higher radiative constant 70 and a larger 
emission wavelength A. From this point of view, J- 
aggregates of squarylium dyes may be suitable candi- 
dates.—'^'^ This type of aggregates, spin-coated on a 
substrate, shows a sharp absorption peak at A « 800 nm 
with HWHM = 20 nm at room temperature and a fast 
(~ 100 fs) optical response^^'^^ combined with a giant 
cubic succeptibility;^ both attributed to the excitonic 
nature of the optical excitations. The monomer decay 
time has been reported to be ^ 100 psj^'^), although no 
information about the quantum yield has been presented. 
If we assume that this time is of radiative nature, the su- 
perradiant constant Tn can be adjusted to values above 
the bistability threshold even for smaller concentration of 
monomers in the film. On the other hand, for larger 70 
also the intensity required for switching increases, which 
is not desired because of the limited photostability of 
most J-aggregates. 

VII. SUMMARY AND CONCLUDING 
REMARKS 

We theoretically studied the optical response of an ul- 
trathin film of oriented J-aggregates, with the goal to ex- 
amine the possibility of bistable behavior of the system. 



The standard Frenkel exciton model was used for a sin- 
gle aggregate: an open linear chain of monomers coupled 
by delocalizing dipole-dipole excitation transfer interac- 
tions, in combination with uncorrelated on-site disorder, 
which tends to localize the exciton states. We consid- 
ered a single aggregate as a meso-ensemble of two-level 
systems, each one composed of an s-like localized one- 
exciton state and its own ground state. The one-to-two 
exciton transitions have been neglected. 

As a tool to describe the transmission properties of the 
film, we employed the optical Maxwell-Bloch equations 
adapted for a thin film. The electric polarization of the 
film was calculated by making use of a joint probabil- 
ity distribution of exciton energies and transition dipole 
moments, properly taking into account the correlation 
properties of these two stochastic variables. The joint 
distribution function was calculated by numerically di- 
agonalizing the Frenkel Hamiltonian and averaging over 
many disorder realizations. 

We derived a novel steady-state equation for the trans- 
mitted signal in terms of the joint distribution func- 
tion, and demonstrated that three-valued solutions to 
this equation exist in a certain domain of the parameter 
space {Tji, cr*), where is the superradiant constant and 
a* is the half-width-at-half-maximum of the absorption 
band. Our approach allowed us to generalize previous 
results3i*32 correctly account for the stochastic nature 
of the exciton energy and transition dipole moment. Us- 
ing the new steady-state equation, we have found that 
the critical value of the so-called "cooperative number" 
rj?./f*,— which governs the occurrence of bistability of 
the film, is higher than obtained before.— Moreover, in 
contrast to Refs. 31 and 32, we have analyzed the switch- 
ing time, which show a dramatic increase for input in- 
tensities close to the switching point. We also found that 
the "cooperative number" Tji/cr* increases with a*, but 
only slightly, varying between 12 and approximately 25 
within a wide range of a* . Estimating the parameters of 
our model for aggregate of polymethin dyes shows that 
these are a promising candidate to measure the effect. 

Finally, we note that also the microcavity arrangement 
of molecular aggregate a^^'^^i^^'^^i^° is of interest for ap- 
plications. During the last decade, organic microcavities 
have received a great deal of attention because of the 
strong coupling of the excitons to cavity photons, leading 
to giant polariton splitting in these devices?^ The recent 
observation of optical bistability in planar inorganic mi- 
crocavities^^ in the strong coupling regime suggests that 
organic microvavities can exhibit a similar behavior. 
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APPENDIX A: ESTIMATES OF QUANTUM 
INTERFERENCE EFFECTS 



dephasing operator, acting as follows: 

(j.|i?*^-*V|0) = . (A2) 



Our approach to the optical dynamics of a single ag- 
gregate was based on the representation of the aggregate 
as a meso-ensemble of two-level systems with their own 
ground states. The model has its origin in the fact that 
the optically dominant exciton states are localized on dif- 
ferent segments and overlap weakly. In reality, however, 
the ground state of an aggregate (all the monomers are 
in their ground states) is common for all (multi- ) exciton 
states. This results in cross-interference of field-induced 
as well as spontaneous transitions. Below, we provide 
estimates of these additional terms and show that in the 
limit of dominant inhomogeneous broadening of the J- 
band, the cross-interference effects can be neglected. In 
our estimates, we will only consider ground state-to-one- 
exciton transitions. 

We start with the equation for the density operator p 



Ho - d£,p 



rad ^ 



(Al) 



where Hq is the exciton Hamiltonian specified in Eq. ([T]) 
and the term —d£ describes the interaction of the aggre- 
gate with the field £ inside the film. R^^^^ represents the 



J 



Here, 71^0 is the (pure) dephasing rate of the exciton tran- 
sition \v) — > |0), excluding radiative decay. These con- 
stants will be considered on a phenomenological basis. 

The operator i?'''*'* governs the exciton radiative relax- 
ation. It is given by (see, e.g., Ref. [63h 

+ p\u){iy'\~2\0){iy\p\i.'){0\] , (A3) 



where j^i, = 70 ( J2n 'Pvn) is the radiative decay rate of 
the population of the uth. state. Furthermore, 7i/,y' = 
If'v 7^ v') describes the quantum interference in the 
radiative relaxation of the I'th and v'th states, resulting 
from the cross-coupling of different decay channels. It 
reflects the fact that a state when decaying, drives 
another state u' and vice versa. If the transition dipoles 
of all the states are parallel, ^j^^i = {'^wly'y' Y^'^ ■ 

Using Eqs. and (fM]) in Eq. ((XT|) . we arrive at the 
following set of equations for the density matrix elements: 



(A4a) 



Pvv' 



-(lEvv' + ^vv')Pvv' ~ ^ X! ^v^"P^"v' ~ ^ X! 7i^"i"Pi/iy" + i {dv^Pl'o - Pvodv'S) , 1^7^ 1^' , (A4b) 



PvQ 



-{ie^ + r^o)Pz.o ~ ^ X! 7i/i^'Pi^'o - i XI Pi^-^'^-i^'^ ~ '^^i^^ (Pi^i^ ~ Poo) , 



(A4c) 



v' ^v 



v'^v 



Poo 



^Pvv = l 



(A4d) 



Here we introduced the notations: Sw' — £v — Sv'i ^w' 
\{lvv + Iv'v') + Ivo + 7i^'o, and r,,o = \lvv + Ivo- 



Equations (|A4[) differ from those used in the two-level 
model, Eq. ([2]), by several terms. Because all the exci- 
ton states have the same ground state, which is reflected 
in the normalization condition (| A4dp . the low- frequency 
coherences are now involved in the aggregate optical dy- 



namics. They are coupled to the populations [Eq. (|A4a[) ] 
as well as to the high-frequency (optical) coherences p^Q 
[Eqs. (|A4b|) and jXIc)) ] via both the cross-coupling of the 
transitions and the field. In addition, the cross-coupling 
also couples the optical coherences pi/o [Eq. (jA4c[) ]. 

In quantum optics of atomic gases, the cross-coupling 
of transitions has been found to be the origin of many 
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fascinating effects, such as narrow resonances, trans- 
parency and gain without ^gopulation inversion (see for 
an overview Refs.[64 65,66.671). as well as bistability at a 
low threshold j^^'^^i^" All these effects, however, require 
specific conditions: equivalent magnitudes of all the 
and the absence of dephasing and inhomogeneous broad- 
ening. Any deviation from these requirements washes out 
those effects. In particular, this happens for J-aggregates; 
below we argue why all the cross-terms in Eqs. (jA4l) can 
be neglected for these systems. 

The contribution of the cross-terms to a given state v 
always comes in the form of a summation over all other 
states v' . The optical dynamics of the system is deter- 



mined by only several dominant states. If N* is the typi- 
cal localization length, there will be N/N* of such states. 
They are spread over the width of the absorption band, 
given by 2a*. Therefore we can estimate the sum un- 
der consideration by {j* /2a*){N/N*) ~ 7o7V/2ct*, where 
7* — ^()N* is the typical radiative rate of optically dom- 
inant states. The materials we consider typically have 
7o 10® s^^ ^ 10~^ cm~^ and 2a* on the order of 
several tens of cm~^. Then, for an aggregate of length 
N = 500 the ratio joN/2a* ^ 0.1. On this basis, we 
neglect all the cross-coupling terms in Eqs. (jA4| and re- 
place the normalization condition (jA4dp for the whole 
aggregate by the one for a single segment, poo + Pi^f = 1- 
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